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Abstract 

We investigate the domain growth and phase separation of hydro dynamically- 
correct binary immiscible fluids of differing viscosity as a function of minority 
phase concentration in both two and three spatial dimensions using dissipative 
particle dynamics. We also examine the behavior of equal-viscosity fluids and 
compare our results to similar lattice-gas simulations in two dimensions. 

I. INTRODUCTION 

Over the last few years, the dissipative particle dynamics (DPD) model of complex fluids 
has received considerable attention. It has matured from its somewhat arbitrary initial 
formulation into a model with a solid theoretical basis. Furthermore, it has been applied 
with considerable success to a large number of computer simulations of complex fluid systems 
such as colloidal suspensions, polymeric fluids, spinodal decomposition of binary immiscible 
fluids, and amphiphilic fluids. DPD also looks promising for simulating multiphase flows 
and flow in porous media, and is now considered a useful technique alongside the other 
complex fluid algorithms: molecular dynamics, lattice-gas automata, and techniques based 
on the lattice-Boltzmann equation. No single technique can yet be applied to all situations, 
and each has different strengths and weaknesses. While molecular dynamics is in principle 
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the most accurate microscopic approach, in practice it is too slow in both its quantum 
(Car-Parinello) and classical forms because of its excessive detail. Discrete mesoscopic 
methods developed from lattice-gas automata have had some success, but they too have 
problems, such as lacking Galilean invariance. The traditional approach of continuum fluid 
dynamics has met with limited success in modeling behavior on the length and time scales 
characteristic of complex fluids. 

In this paper we investigate phase separation from both symmetric (critical) and 
off-critical quenches in binary immiscible fluids of generally differing viscosity using the 
hydrodynamically-correct DPD model. Our motivation is to probe and extend knowl- 
edge of the behavior of differing-viscosity fluids, of application, for example, to the action 
of detergents and the extraction of oil from reservoir rocks. Phase separation in equal- 
viscosity binary immiscible fluid systems has been simulated using a variety of techniques, 
including DPD [0-0] ; molecular dynamics @-||; Monte Carlo IjSjpTOI; cell dynamical systems 
without hydrodynamics |TT| and with Oseen tensor hydrodynamics [O; time- dependent 



Ginzburg-Landau models with |13| - |17| and without hydrodynamics [fL7j-|19l; lattice-gas 
automata p0^p7|; and lattice-Boltzmann techniques ]f28|-|3l|j. Spinodal decomposition of 
differing-viscosity immiscible fluids has previously been simulated in two dimensions by a 
time-dependent Ginzburg-Landau model without hydrodynamics [J35| . We discuss the ef- 
fect of the proportion of each phase on the scaling behavior in both two and three spatial 
dimensions. We also examine the behavior of equal-viscosity fluids, comparing our results 



to similar lattice-gas simulations in two dimensions [27 



After describing our fluid model in the following section, we discuss the expected temporal 
development of the characteristic size of the separating domains in Section [TTI|. We then 
describe our method for calculating the characteristic domain size and its rate of growth 
in Section [TV]. This is followed in Sections [V] and [VT| by information on the simulations 
performed and a discussion of the results, and by some conclusions in Section [VI]] . Finally, 
as an Appendix to this paper, we make a few comments on the high-performance computing 
aspects of this work. 



II. THE FLUID MODEL 

In 1992, Hoogerbrugge and Koelman proposed dissipative particle dynamics |56] as a 
novel particulate model for the simulation of complex fluid behavior. DPD was developed in 
an attempt to capture the best aspects of molecular dynamics and lattice-gas automata. It 
avoids the lattice-based problems of lattice-gas automata, yet maintains an elegant simplicity 
and larger scale that keeps the model much faster than molecular dynamics. This simplicity 
also makes DPD highly extensible, such as for including the interactions of complex molecules 
or modeling flow in an arbitrary number of spatial dimensions. The key features of the model 
are that the fluid is grouped into packets, termed "particles" , and that mass and momentum 
are conserved but energy is not. Particles are normally interpreted as representing a coarse- 
graining of the fluid, so that each particle contains many molecules |36|-p9|. Since the 



intrinsic time scale of DPD represents the correlated motion of mesoscopic packets of atoms 
or molecules, it is typically orders of magnitude larger than the time scale of molecular 
dynamics PD |. Particle positions and momenta are real variables, and are not restricted to 
a grid. 
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Espanol and Warren's analysis |37] showed that the original DPD model does not satisfy 



detailed balance, so the equilibrium states (if they exist) cannot be simply characterized. 
Detailed balance is the condition equating the rates of forward and backward transition 
probabilities in a dynamical system, and is a sufficient (but not necessary) condition guar- 
anteeing that the system has a (Gibbsian) equilibrium state fH|,fiS| . Espanol and Warren 
formulated a Fokker-Planck equation and equivalent set of stochastic differential equations 
which lead to an analogous continuous-time model, 

d Pi = J2 Fydt = [ F ?j dt + F S dt + v u' nv U 

j¥« i^i (1) 

<ix 7 - = M^dt. 
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In these equations, Pi, x i; and rrii denote the momentum, position, and mass of particle i. Fg 
is a conservative force acting between particles i and j, while Fg and Fg are the dissipative 
and random forces. dWij = dWji are independent increments of a Wiener process. By Ito 
calculus 

dWijdWki = (SikSji + SaSjk) dt, (2) 
so dW^ is an infinitesimal of order ~ and dW^ can be written Oijy/dt, where % = 9ji is a 



random variable with zero mean and unit variance |J0]| . Detailed balance is satisfied by this 
continuous-time version of DPD with an appropriate choice for the form of the forces |43| . 
and so equilibrium states are guaranteed to exist and be Gibbsian. To ensure that the 
associated fluctuation-dissipation theorem holds, Espanol and Warren suggested the forces 
assume the following forms [B7|]: 

Fg = aWijGij, (3) 



F 5 = -7<4(e<j ■ v ii)e*i, (4) 

and 

Fg = aUijBij, (5) 

where = (Pi/wj) — (pj/rrij) is the difference in velocities of particles j and i, is the unit 
vector pointing from particle j to particle i, and Uij is a weighting function depending only 
on the distance separating particles i and j. The constants a, 7, and a are chosen to reflect 
the relative importance of the conservative, dissipative (viscous), and random components 
in the fluid of interest. As a consequence of detailed balance and the fluctuation-dissipation 
theorem, 7 and a are related to Boltzmann's constant kg and the equilibrium temperature 
Tby 

2 

— = 2k B T. (6) 
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In order to remain as close as possible to the original DPD model, Espanol and Warren 
chose the friction weight function to be 
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(7) 



within the constant cutoff length r c > and zero otherwise, where is the distance between 
particles i and j. Adding Eqs. (|3|)-(|5]) together, the total force is 



F- 



a - -yuij (eij ■ v i3 -J + 



Uifiiy (8) 



Immiscible fluid mixtures exist because individual molecules attract similar and repel 
dissimilar molecules. The most common example of such behavior arises in mixtures of oil 
and water below a critical temperature. The nonpolar, hydrophobic molecules of oil attract 
one another through short-range van der Waals forces, while the polar water molecules have 
complex, long-range hydrophilic attractions which are dominated by electrostatic interac- 
tions, including hydrogen bonds. At the atomistic level employed in molecular dynamics, 
such interactions demand a detailed treatment. However, to obtain mesoscopic and macro- 
scopic level descriptions using DPD, the microscopic model can be drastically simplified. 

In order to model immiscible fluids, the simplest modification to the one-component 
dissipative particle dynamics algorithm is to introduce a new variable, called the "color" by 



analogy with Rothman-Keller p0| , |36| . For example, we could choose red to represent oil 
and blue to represent water. When two particles of different color interact we increase the 
conservative force, thereby increasing the repulsion. That is, 

J «o if particles i and j are the same color , . 

%d \ ol\ if particles % and j are different colors 

where «o arid «i are constants with < «o < a i- As a consequence of mass and momentum 
conservation, the Navier-Stokes equation is obeyed within a single-phase DPD fluid and 



within regions of homogeneity of each of the two binary immiscible fluid phases p6| , p7| , |43| , P4 . 
Likewise, detailed balance is also preserved, at least in the limit of continuous time p7l , |43l , |4"4" . 

The immiscible fluids described above are identical to each other. However, it is often 
the case that the two fluids in a mixture will differ physically. For example, oil and water 
typically have different viscosities. To model binary fluids with differing viscosity we again 
adopt the simplest approach: labeling one of the two phases (colors) as more viscous. When 
two particles of the same viscous color interact, the dissipative (viscous) force is increased; 
in order to keep the temperature constant, we must correspondingly decrease the random 
force according to Eq. (||), i.e., 

J 7o if either particle i or j is not the viscous color , , 

^ ^ \ 7x if both particles i and j are the viscous color. 



where 



, do if either particle i or j is not the viscous color . . 

' ai if both particles i and j are the viscous color, 



U l = °1 = 2k B T. (12) 
7o 7i 
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Our previous study of finite difference methods applicable to simulations with non- 
conservative forces Q indicated that the finite-difference algorithm suggested by Groot 
and Warren [HTJ] is a good choice for DPD. Their method is a variation on the familiar 



velocity-Verlet algorithm, adding a momentum estimate before the force evaluation: 



x 4 (t + At) =Xi(t) + 



At 

rrii 
At 



Pi (t) + 



Pi(t + At)=p i (t) + —i i (t) 
i i (t + At) = J EF ij (t + At) 

Pi (t + At) = Pi (t) + %(t) +fi(t + At)} 



(13) 



where At is the time step size and fj(t) is the force acting on particle i at time t. The DPD 
units of length, mass, and time are specific to the particular set of model parameters, and 
the exact relationship between these parameters and a real fluid in a particular situation 
can be determined by considering the dimensionless groups relevant to the motion of that 
system, such as the Mach, Reynolds, and Weber numbers. 

Many further modifications to the model have been suggested and implemented by others. 
Some of the most interesting include energy conservation, colloidal particles, and polymers. 
In very recent work, it has been shown that it is possible to derive a modified version of 
DPD directly from the underlying molecular dynamics 



III. TEMPORAL BEHAVIOR OF THE CHARACTERISTIC DOMAIN SIZE 



A central quantity in the study of growth kinetics is the characteristic domain size R(t). 
For binary systems in the regime of sharp domain walls, this is usually thought to follow 
algebraic growth laws of the form 



R(t) ~ t 13 . 



(14) 



For symmetric quenches without hydrodynamic interactions, dynamical scaling theory and 
experiment j|7| indicate that the scaling exponent /3 = ~. If flow effects are relevant, 
typically 



in two dimensions, and 



P 



for R<t:R h 
for i? > R h 



for < R d 

for R d < R < R h 

for i? > R h 



(diffusive) 

(inertial hydrodynamic) 



(diffusive) 

(viscous hydrodynamic) 
(inertial hydrodynamic), 



(15) 



(16) 



in three dimensions, where Rh = i] 2 /(pn) is the hydrodynamic length and R d = y/rjD is 
the diffusive length, expressed in terms of the absolute (dynamic) viscosity r/, density p, 
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surface tension coefficient k, and diffusion coefficient D. These scaling laws follow directly 
from dimensional analysis of the macroscopic fluid dynamics equations (so-called model H, 
or Cahn-Hilliard coupled to Navier-Stokes hydrodynamics) in the appropriate regimes [37 



When R <C Rd in three dimensions, diffusive effects dominate and we expect the domains 



to form via the Lifshitz-Slyozov- Wagner (LSW) evaporation-condensation mechanism P? 
When R 3> Rh in both two and three dimensions, we expect the growth mechanism to 
be surface tension driven by hydrodynamic flow, balanced by inertial effects [38|. If Rd <C 
R <C Rh in three dimensions, we expect viscous hydrodynamic effects to dominate, as 
predicted by Siggia [33 1; in this regime the surface tension drives the transport of the fluid 
along the interface, which is possible only if both phases are continuous P7|j49[] . This third 
growth regime cannot occur in two dimensions, and so we expect to see diffusive growth for 
R -C Rh [f38 |. Due to our choice of model parameters and the small size of our simulations, 
we do not expect to probe the viscous or inertial hydrodynamic regimes. 

Simulations in two spatial dimensions are especially useful to emphasize the importance of 
correct hydrodynamics: simulations which do not conserve momentum typically give (3 = ~ 
for the diffusive regime (R <C Rh) PH3|P^p6|- Simulation methods with correct hydrody- 
namic interactions typically give the expected result of (3 = | (see e.g., Refs. [pl| |3|j6| JT|,p~0|,p~6f| ) . 
It is worth noting that momentum conservation is not thought necessary to model spinodal 
decomposition in binary metallic alloys, since phase separation occurs by the migration of 
atoms to neighboring vacancies on the crystalline lattice [50|. Simulations based on lattice- 
Boltzmann techniques also typically display /3 = |. These observations are supported by a 
renormalization group approach which has shown that thermal fluctuations cause Brownian 
motion-driven coalescence and play a crucial role in causing (3 to assume the value of ^ p? 



29 , most do not. 



although some lattice-Boltzmann techniques include these fluctuations 

There is some confusion in the literature over which scaling exponents should be observed 
for off-critical quenches, and whether or not the algebraic scaling law [Eq. (0)] should in 
fact be obeyed for any off-critical binary immiscible fluid. Several authors have reported the 
coexistence of multiple length scales 



It is likewise uncertain as to what behavior 

Certainly, one 



54,51 



should be observed from binary immiscible fluids of differing viscosity f35 



growth mechanism we expect for off-critical quenches of both equal and differing viscosity 
is the LSW evaporation-condensation mechanism, for which (3 = | |47] . 

It should also be noted that there are still some experimental and theoretical challenges in 
unraveling the behavior of systems in which both the order parameter and the momentum 
are locally conserved W7\. Experimentally, for example, it is difficult if not impossible 
to study two-dimensional fluid systems. As far as numerical studies are concerned, it is 
important to recognize that three-dimensional simulations are particularly demanding on 
all the aforementioned techniques, and so definitive results are harder to come by than in 
two dimensions. 

Indeed, Cates's group in Edinburgh has recently reported somewhat conflicting results 
from three-dimensional studies of binary immiscible fluid separation for equal viscosity flu- 
ids using dissipative particle dynamics and lattice-Boltzmann methods While the 
former suggests the persistence of non-universal length scales, hypothesized to be due to 
the intrusion of a "molecular" or discretization length scale, this is not supported by the 
latter, from which finite-size effects were more rigorously excluded. Note, however, that 
whereas the lattice-Boltzmann scheme was based on a Landau free-energy approach, and is 
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essentially no more than a finite-difference solution of the continuum model H equations, 
the macroscopic equations to which the spinodal DPD scheme corresponds have not yet 
been derived. This makes it unclear whether the two systems being simulated are really 
one and the same. Both of these studies emphasize the importance of observing dynamical 
scaling over decades of time before drawing conclusions as to the true nature of the scaling 
regime. One notable result of their work is the clarification of the relative time scales over 
which computer simulation techniques typically operate. For the simulations they discuss, 
molecular dynamics time scales are on the order of 10 2 in dimensionless units, lattice-gas 
automata and Langevin dynamics probe time scales on the order of 10 2 through 10 4 , DPD 
typically 10 3 to 10 5 , and methods based on the lattice-Boltzmann equation 10 1 through 10 8 . 

Grant and Elder have recently argued |53| that (3 < 5 in any asymptotic scaling regime 
because the Reynolds number Re = pRV/r] (where V is a characteristic velocity) cannot 
diverge, and in fact must remain less than a critical value Re cr to avoid the onset of turbu- 
lence |j3 1 and possible turbulent remixing of the fluids. The conclusion they draw is that 
the P — I scaling regime must be transient |53[. However, their analysis neglects men- 
tion of the relative strength of the interface, quantified for example by the Weber number 
Nwe = pRV 2 /n. If iVwe is small at the onset of turbulence we would expect turbulent 
remixing to be delayed, or perhaps even postponed indefinitely allowing Re to diverge. In 
any case, the separation dynamics would likely be affected by Re > Re CT ; for example, the 



fluid viscosity in a turbulent region is roughly proportional to Re [54 



IV. ESTIMATING THE CHARACTERISTIC DOMAIN SIZE 

In order to characterize the phase separation kinetics within a binary immiscible fluid, 
we need a practical tool to measure the characteristic domain size corresponding to the 
state of the system at a given point in time. The use of the static structure function for 
this purpose is widespread. However, we have noticed bizarre "early time" behavior when 
using the static structure function to characterize simulations of highly off-critical quenches. 
Specifically, the characteristic domain size would sometimes change suddenly by an order of 
magnitude. This abrupt change in behavior did not correspond to anything observable in 
the time evolution of the positions of the DPD particles. Such anomalous behavior is likely 
due to the large fluctuations prevalent in the static structure function for small simulations 
of highly-mixed binary fluids, for which the characteristic domain size is very small. 

The radial distribution function g(r) is a well-established tool for the analysis of single- 
phase fluids |55| , |56[] , and indicates the likelihood of finding two particles separated by a 
distance r. For binary fluids we can also use the same-phase and differing-phase distribution 
functions ||. The same-phase distribution function <?oo( r ) describes the likelihood of finding 
two particles of the same phase separated by a distance r, and the differing-phase distribution 
function g Q1 (r) describes the likelihood of finding two particles of differing phase separated 
by a distance r. From the peaks of the differing-phase radial distribution function, we can 
estimate the characteristic separation of particles of differing phase (i.e., the characteristic 
domain size). Consequently, we decided to calculate the distribution between particles of 
different phase (color), 
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^ 0l(r) " piV0[l-0]V(r,Ar)' (17) 

where m = Vi is the mass of each particle (assumed all identical), n(r, Ar) is the number 
of particle pairs of differing phase with separation between r — Ar/2 and r + Ar/2, N is 
the total number of particles, G (0, 1) is the fraction of particles of one phase (the more 
viscous phase if the two phases are of different viscosity), and V(r, Ar) is the volume of the 
spherical shell of radius r and thickness Ar (from r — Ar/2 to r + Ar/2). It is apparent 
from Eq. ( |17D that (?oi( r ) can only be calculated for r < L/2 in a simulation with a periodic 
box size of L; this should not be interpreted as an undesirable limitation as finite size effects 
are normally significant for R > L/2 PJ5"2"| . A value of Ar = L/2000 gives a reasonable 
compromise between noise and resolution for the size of simulation we discuss in this paper. 

The principal difficulty lies in analyzing the results automatically, as we need to calculate 
9oi{ r ) a t many time steps within each simulation in order to display the domain growth over 
time. We chose to calculate the median of that portion of the smoothed goi ( r ) curve extend- 
ing above a threshold value 1 + 3s, where s is the standard deviation of the smoothed curve 
when it has effectively a constant unit value, estimated from the last tenth of the smoothed 



goi(r) curve at t = 1. We used a fourth order Savitzky-Golay smoothing filter |]57| , p8| to 
smooth g i (r) over r at each point in time with a symmetric window size chosen to reduce 
the noise while leaving significant features untouched (41 points in two dimensions and 101 
points in three). We chose the cutoff threshold of 1 + 3s so as not to bias the median by the 
size of the periodic simulation box; likewise, we chose the median in preference to the mean 
(first moment) because the median is less biased by outliers, such as those 2.5% of noisy 
points which effectively have unit value but extend above the threshold. Both the median 
and mean give poor estimates of the characteristic domain size as it approaches L/2, the 
limit of calculable goi (r). When this situation is detected, we use the global maximum of 
goi( r ) to estimate the domain size instead; this is a continuous transition for symmetrical 
peaks. 

Once we have calculated the characteristic domain size for a series of goi ( r ) curves taken 
at different times from a single simulation, we begin a search for linear sections in the plot 
of log 10 R vs. log 10 1. This has also been automated, using analysis of variance to decide 
whether a given section of the plot is linear or cubic. We used the analytic expressions for 
the least-squares fits with moments taken about the means to minimize the effect of roundoff 
error. We keep only the longest linear sections, subject to there not being any significant 
gaps in the coverage of the log-log plot. An ensemble average of a number of simulations 
yields a plot of the scaling exponent (3 vs. log 10 1, in which long horizontal (zero gradient) 
sections represent algebraic growth. This procedure allows more accurate determination 
of the scaling exponent than a visual comparison of log 10 R vs. log 10 1 to a straight line of 
a particular slope, and provides a statistically valid determination of whether or not the 
growth is truly algebraic over a particular time period. Finite size effects for R m L/2 
normally result in non-algebraic growth or unusually low scaling exponents; either case is 
easily detected by our method. A further advantage is the automation we have described, 
permitting large ensemble averages with minimal effort. 

We should note that estimating the characteristic domain size using the median can 
occasionally lead to discontinuities. If these discontinuities are large enough they will cause 
a break between linear sections. Small discontinuities may be spanned by a single line, which 



S 



would affect the slope (ft) and be detrimental to the results. However, large discontinuities 
will not be spanned, and so will not affect the slope. 

The results we obtained using these techniques are considerably better than those ob- 
tained with the static structure function, especially for highly-mixed fluids where the pair 
distribution function is not drowned by fluctuations to the extent that the static structure 
function is. Furthermore, the pair distribution function is more intuitive to analyze than 
the static structure function, which aids our interpretation of the results. 



V. SIMULATION RESULTS 

For the simulations of fluids with differing viscosity, we chose one phase to be an order of 
magnitude more viscous than the other, i.e., 71 = IO70. Before simulating this new complex 
fluid, it is advisable to verify that increasing the parameter 7, while keeping the temperature 
constant, does indeed increase the viscosity. For the lower- viscosity fluid we chose the model 
parameters shown in Table |. The absolute (dynamic) viscosity of a homogeneous DPD fluid 
can be estimated theoretically from the continuous-time viscosity, ignoring the effect of the 
conservative forces (i.e., a=0) |59|j60[| ; the continuous-time viscosity of this lower- viscosity 



fluid is 7] = 2.8. 

In order to verify this estimate, we performed a series of simulations of steady shear 
using Lees-Edwards periodic boundary conditions. We performed a total of 63 simulations 
with a time step of At = 0.1 (in our DPD time units), each from a different random initial 
configuration and each allowed to settle to steady-state shear before measurement began. 
We studied systems of both 1600 and 6400 particles, six simulations at each of nine different 
shear rates for the former and three simulations at each of three distinct shear rates for the 
latter. As the results from the larger simulations gave a mean viscosity nearly identical to 
that of the smaller ones, we can conclude that finite size effects do not bias the viscosity of 
the smaller, faster simulations. We calculated the velocity profile for each set of parameters, 
and found it to be statistically indistinguishable from linear in every case. 

Analyzing these simulations led to a conclusion of 77 = 1.94 ±0.01. Others have also 
found discrepancies between theory and simulation, particularly regarding the kinematic 
contribution to viscosity |5{|, so the difference between the simulated viscosity and the 
continuous-time viscosity is not entirely surprising. Since molecular dynamics simulations 
containing only conservative forces give a finite viscosity, we would be surprised if the theo- 
retical estimate did not differ from our calculations. 

In order to measure the viscosity of the more viscous fluid, we set up a series of 1600 
particle simulations of a homogeneous DPD fluid as described above, using a total of 30 
simulations. We varied the shear rate in the range that gave the best results in the previous 
simulations. The model parameters differed in that 7 was a factor of ten larger while a was 
a factor of \/T0 larger to keep the temperature constant [see Eq. (|i~2"|)l; these parameters are 
shown in Table [n[ We decreased the time step to At = 0.01 due to the increased magnitude 
of the dissipative and stochastic forces. Analyzing these results led to the conclusion that 
T) = 17.2 ± 0.3, which confirms the increase in viscosity. This increase is close to a factor of 
ten, so that for similar model parameters it is reasonable to conclude that 7 is approximately 
proportional to viscosity. 
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It is possible to measure the surface tension of a binary immiscible fluid by integrating the 
pressure tensor across a flat interface, or by verifying Laplace's law for a series of equilibrium 
bubbles of varying radii ||55|| . Laplace's law was verified for a DPD binary immiscible fluid in 
our previous simulations |]-f3|, and so surface tension measurements were omitted from the 
present study. These calculations would, however, allow identification of the unit of time for 
comparison with other simulation techniques (see e.g., Refs. [|5l,|52| ) . Theoretical estimates 
for surface tension are also available, but are of similar accuracy to the viscosity estimate 
above. 

For the simulation of fluids of differing viscosity, we must also choose the relatively small 
time step size of At = 0.01 in order to ensure the stability of the algorithm as a result of 
the increased magnitude of the dissipative forces. This has the consequence of making these 
simulations computationally much more expensive than equal-viscosity simulations. In this 
context, it is worth commenting on the virtues of using DPD to perform simulations of 
differing-viscosity fluids compared with other models. Sappelt and Jackie use an approach 
based on the Cahn-Hilliard equation (so-called model B without noise) with a concentration- 
dependant mobility |35| so do not include hydrodynamics, unlike our approach with DPD. 
The other mesoscale techniques which could be used to model these fluids include lattice- 
gas automata and methods based on the lattice-Boltzmann equation. As with DPD, there 
is a high computational price to the lattice-gas approach, which requires adjustment to 
the collisional outcomes of the look-up tables. Lattice-Boltzmann (more correctly, lattice- 
BGK) models can be used, based for example on the Swift-Osborn-Yeomans free energy 
functional approach ||61|| , but there is a similar computational price, although with the 
additional complication of poorly understood numerical instabilities. However, DPD offers 
the simplest algorithmic implementation which is thermodynamically consistent. 

Our differing-viscosity simulations used the model parameters shown in Table |TTT]. Each 
simulation had 6400 particles, and ran for 50,000 time steps. We chose to use only 6400 
particles in our simulations so that individual simulations would be quick enough to be 
run multiple times, allowing us to consider the effect of changing the phase fraction and 
viscosity, and to increase the confidence in our results and calculate accurate error estimates 
with ensemble averaging. Simulations with more particles would have given better resolution 
of small scale features (relative to the system size) and postponed the finite size effects, at 
the price of increased computational demands. Our resources would only have allowed a 
few simulations of the size and length used by Jury et al. || (1,000,000 particles), making 
studies of the sort we describe in this paper impossible. 

We calculated the pair distribution function at 64 logarithmically-spaced points in time 
during the course of each simulation, starting from time t = 1 and finishing at t = 500 
(in our DPD time units). In both two and three dimensions, we performed ten simulations 
at each of nine different fractions of the viscous phase, ranging from <fi = 0.1 (10% viscous 
phase) to <p — 0.9 (90% viscous phase). We show the time evolution of a single simulation 
at each value of the viscous fraction in Figs. [[]-§] for both two and three dimensions. We 
represent the positions of the DPD particles in two dimensions by a gray scale map, in 
which the particle positions are weighted by Eq. (|7D and assigned an intensity of gray based 
upon the proportion of each phase in this localized-average. Figures |3] and |] display the 
three-dimensional surface of the interface between the two immiscible phases, as defined by 
there being equal proportions of each fluid in the localized-average. In these four figures the 
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gradual development of domains can be seen. We examined the results of the domain size 
analysis for each simulation in plots of log 10 R vs. log 10 1, and further examined each ensemble 
average in plots of (3 vs. log 10 1; the latter yielded the mean values of (3 and corresponding 
68% confidence intervals for significant sections of algebraic growth. We show these results 
for two and three dimensions in Tables |V| and [V], and Figs. [5] and || The range of log 10 1 
we give in the tables should be taken as a rough guide, since there is often a high degree 
of correlation between the start and end times of a particular linear section and its growth 
exponent within the broad category of late time. We should emphasize that we are using 
the term "late time" as a relative category in this paper, to distinguish these results from 
those obtained at very early times. Due to the model parameters and small size of these 
simulations, we do not expect to probe the viscous or inertial hydrodynamic regimes (see 
Eqs. dID-(|§), and compare with Refs. |Tj-[§). 

We also constructed a series of equal-viscosity simulations in both two and three dimen- 
sions. These had At = 0.1, 6400 particles, and stopped at t = 1000. The model parameters 
were the same as in the differing-viscosity simulations, with the exception that both phases 
were identical to the less- viscous phase of the differing-viscosity simulations (i.e., 7 = 70 
and a = 00). We calculated the pair distribution function at 64 different points from t — 1 
to t — 1000, with ten simulations at each of eight different minority phase fractions from 
<p = 0.05 to cf) = 0.5 in two dimensions, and ten simulations at each of five different fractions 
from cf) = 0.1 to cf) = 0.5 in three dimensions. We show the time evolution of a single simu- 
lation at each value of the minority fraction in Figs. |7] and |8| for two and three dimensions 
respectively. In these two figures the gradual development of domains can be seen, and 
is greatly slowed for small minority fraction. We show the scaling exponents for two and 
three dimensional equal- viscosity fluids in Tables VI and VTJ, and Figs. |^ and 0, where we 
mirror the scaling exponents about = 0.5 to show the full range of minority phase fraction. 
Comparison with Figs. [7] and || qualitatively confirms the same behavior. 



VI. DISCUSSION 

A feature common to all these simulations is the lack of scale-invariance at very early 
times, until approximately t — 10 to t — 50 (in our DPD time units) depending on the 
exact composition of the fluid. This is apparent in the goi ( r ) curves as multiple peaks 
of similar magnitude, but could not be seen in the frequency domain shown by the static 
structure function. Figure [TTJ shows an example of the differing-phase radial distribution 
function for a three-dimensional equal- viscosity simulation (<f) = 0.5) at t = 13.9; the crosses 
represent the actual data and the curve represents the smoothed data. This is one of the 
simulations we show in Fig. |8|. The multiple peaks typically evolve by changing their height 
and relative weight, but not their position. This is noticeably different from our "late time" 
behavior of the distribution function, where a single peak gradually advances its position 
while broadening but retaining its height and general shape. It is because of this "early time" 
behavior that it was decided not to use solely the global maximum of goi (r) to estimate the 
characteristic domain size. 

At these "early times" in all of the simulations, we observed algebraic growth with a 
very small exponent, roughly 0.126 ± 0.003 in two dimensional equal- viscosity simulations, 
0.062 ± 0.009 in two dimensional differing-viscosity simulations (except <p = 0.9 (90% vis- 
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cous phase), for which (3 = —0.24 ± 0.07), and 0.0246 ± 0.0008 for simulations in three 
dimensions. This corresponds to the region of breakdown of scale invariance we described 
at the beginning of this section. These "early time" exponents are unaffected by viscous or 
minority phase fraction in all of the simulations, except the differing-viscosity simulations 
in two dimensions. Here f3 decreases with increasing (viscous phase fraction), probably 
due to the domain-growth arresting effect of increasing viscosity. We observed a remarkable 
(though short) regime of (3 = —0.24 ± 0.07 for = 0.9 (90% viscous phase), for which we 
have no adequate explanation. Others normally discard similar early time regimes without 
comment [^,|T^,^,^,^,^,^] or as an "early stage" or "transient" regime PP^. However, 
there is growing evidence for the coexistence of multiple domain sizes and hence a breakdown 
in universality, at least in certain phase-ordering domains [p^ ,j^, |5l | l . 

For "late time" domain growth in the two-dimensional differing-viscosity simulations 
(see Fig. ^|), we observed a fairly constant value of (3 for = 0.2 (20% viscous phase) 
through = 0.6 (60% viscous phase), decreasing both for = 0.1 and very slightly for 
= 0.7 and = 0.8, then decreasing sharply at = 0.9. This asymmetry is consistent 
with the variation of the "early time" exponent, in that an increasingly viscous fluid is 
expected to develop domains more gradually. At increasingly rarefied fractions (0 = 0.1 
and = 0.9), domain growth is retarded by the increased isolation of the droplets. The 
"late time" growth exponent throughout is effectively |, which suggests that the presence 
of fluids of differing viscosity interferes with the normal (3 = | growth mechanism in two 
dimensions. The growth exponent of | is expected from the LSW evaporation-condensation 



mechanism J47J. This is in some ways analogous to the effect obtained by deliberately 
breaking momentum conservation in symmetric quenches, as described previously 01 • 
Our domains are considerably less circular at all viscous fractions than those observed by 
Sappelt and Jackie (compare Figs. [I] and Q with Ref. J35J). This is likely due to the lack of 
hydrodynamic interactions in their model and to the greater difference in viscosity between 
their two phases. As in their simulations, our two-phase structure for fluids of differing 
viscosity is not very different from the structure for fluids of equal viscosity. Moreover, our 
simulations do not reveal any new insights regarding interfacial structure. 

In three dimensions, the "late time" domain growth of differing-viscosity fluids displays 
nearly the opposite behavior, with the scaling exponent increasing as the viscous fraction 
reaches its extremes. This could be explained by the increased fluid mobility in simulations 
with an extra spatial dimension, as the majority phase is completely connected and so 
the domain growth could occur according to the (3 = | mechanism, which is surface tension 
driven by hydrodynamic flow, balanced by inertial effects. This is qualitatively substantiated 
by inspection of Figs. [3]and|], where a larger degree of connectivity of the majority phase can 
be observed at the extremes of viscous fraction than for « 0.5. However, a more obvious 
mechanism for the domain growth would be the (3 = | LSW evaporation-condensation 
mechanism [ |7|] ; it may be a combination of these two mechanisms that leads to our observed 
|</3<! growth. A slight asymmetry in the "late time" growth exponent is also evident, 
with domain growth proceeding more slowly with increasingly viscous fluids of non-extreme 
viscous fraction. 

For equal-viscosity fluids in two dimensions we observed the expected /3 = | for sym- 
metric quenches (0 = 0.5) |j]-|3[]. As we reduced the minority phase fraction, we observed 
a steady decrease in the scaling exponent until j3 = | is reached at the extremes. This 
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confirms the results of other workers |27] that increasingly off-critical quenches retard the 



domain growth, while providing support for the observed slowdown of growth at the extremes 
of viscous fraction in differing-viscosity fluids in two dimensions, which we commented on 
above. 

In three dimensions, the domain growth of equal-viscosity fluids appears largely unaf- 
fected by varying the minority phase fraction. Although there may be some increase in 
j3 at the extremes of (p (as seen in the differing-viscosity fluid in three dimensions), this is 
difficult to confirm definitely because of the large variation in rate of growth observed for the 
simulations with = 0.1, and hence correspondingly large confidence interval. The scaling 
exponent throughout is close to (3 — |. Whereas Jury et al. || were intending to probe the 
viscous or inertial hydrodynamic regimes with their DPD simulations of equal- viscosity fluids 
in three dimensions, we aimed only to probe length scales below Rd and Rh [see Eqs. (^4|)- 
(fyf)]. As such, our results are fully consistent with theirs. Our exclusion of finite size effects 
is more rigorous than theirs, and although not as extreme as that advocated by Kendon et 
al. [|5^] our method gives statistical confidence that these domains are scaling algebraically. 
Both Jury et al. || and Kendon et al. were able to cover the time domain more fully in 
three-dimensional equal- viscosity symmetric quenches only at the cost of performing a large 
number of computationally very intensive and very expensive simulations. 



VII. CONCLUSIONS 

In this paper, we have described simulations of the domain growth and phase separation 
of hydrodynamically-correct binary immiscible fluids of differing and equal viscosity as a 
function of minority phase concentration in both two and three spatial dimensions. Due to 
our choice of model parameters and the small size of our simulations, we did not expect 
to probe the viscous or inertial hydrodynamic regimes. In three dimensions, we found that 
the characteristic domain size scales as t 1 ^ 3 for simulations of differing and equal-viscosity 
fluids developing from symmetric and slightly off-critical quenches. For highly off-critical 
quenches we observe an increase in the scaling exponent. In two dimensions, we also observe 
t 1//3 in simulations of differing-viscosity fluids developing from symmetric and slightly off- 
critical quenches, although we observe a decrease in the scaling exponent for highly off- 
critical quenches. In equal-viscosity fluids in two dimensions, we observe t l l 2 for symmetric 
quenches and a roughly linear decrease to t 1//3 for highly off-critical quenches; these results 
are in agreement with similar lattice-gas simulations in two dimensions | |27| |. 

Obtaining meaningful results for ensemble averages of highly off-critical binary immisci- 
ble fluids was only made feasible by our automation of the calculation of the characteristic 
domain size by the pair correlation function. It also made possible the identification of a 
regime of breakdown of scale invariance at very early times, which was not noticeable in our 
original analysis using the static structure function. Further simulations aimed to probe the 
viscous and inertial hydrodynamic regimes [see Eqs. (|T^)-(|l6|)1 would be a useful addition to 
this work, as would simulations aimed to cover longer periods of time; however, both would 
require substantially increased computational work. 
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APPENDIX: HIGH-PERFORMANCE COMPUTING 

In this appendix we provide a few comments regarding the running of DPD simulations on 
high performance computers. We usually had easy access to single-processor workstations, 
with a large variety of types and speeds of processors. Multi-processor machines allowing 
parallel execution of simulations are much less common and are more difficult to obtain 
access to, although they have become more common during this research project. We used 
both the Cray T3D of the Edinburgh Parallel Computing Centre (EPCC) and the Hitachi 
SR2201 of the Cambridge High Performance Computing Facility (HPCF) for computing the 
results described in this paper; the former consisted of 512 processor nodes and the latter 
consists of 256 nodes. 

The implementation of the dissipative particle dynamics algorithm is very similar to that 
of conventional molecular dynamics algorithms |p5]| . For example, we divide the periodic 
spatial domain (the simulation cell) into a regular array of equally sized link-cells, such that 
each side of the rectangular domain has an integer number of cells and each cell is at least 
r c across. Each link-cell consists of a dynamically allocated array of particles, and pointers 
to the neighboring cells. Individual particles consist of the position-momentum vector pair 
and a color index. 

For each time step we iterate through the particles in each link-cell, calculating the 
force acting on each particle as it interacts with the particles in the same and neighboring 
link-cells. Since the DPD force acts between pairs of particles, we must ignore half of 
the neighboring cells to avoid duplication. When considering a different particle pair, we 
compare the square of the separation distance with skipping to the next particles if the 
pair is out of range. We then compute the new position and velocity as determined by the 



finite-difference algorithm (see Eq. (13) and Ref. [[|). 

We may write the complete state of the system to file, and we can perform other cal- 
culations thereafter, for example to determine the temperature and pressure of the system. 
We used the freely-available Gnu-make utility to dictate the compilation process, since the 
decision structures it contains make it simple to write programs portable to a large range 
of architectures. We created a comprehensive, automated test suite to make it easy to 
verify that optimizations of the calculations did not accidentally change the results of the 
computations. 
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Given constant r c and number density n — p/m, the DPD algorithm scales linearly (in 
both computation time and memory size) with increasing number of particles (N), and is 
limited by computation time on all but the smallest machines. The main simulations we 
performed for this paper consisted of 6400 particles, and it is on the parallel and serial perfor- 
mance of this size of simulation that we will make most of the following comments. Details 
of the performance of this size of simulation in two dimensions are shown in Tables [VIII 
and [TX|. These tables give the elapsed time per node in seconds and relative parallel effi- 
ciency for the first 1000 time steps, including data for a variety of computers and partition 
sizes. These data are for code compiled with the highest level of optimization, including 
some small reductions in floating-point accuracy. Table |V111| describes the computers used 
to calculate the results in this paper, while Table |IX| describes the computers to which we 
have recently been allowed access, such as the Computer Services for Academic Research 
(CSAR) Cray T3E and SGI Origin2000 in Manchester. 

A typical simulation of 50,000 time steps takes 2.5 hours on a 350 MHz Intel Pentium II 
PC, the fastest single-processor machine to which we had common access. This same sim- 
ulation would take 2.9 hours on a 16-node partition of the T3D at the EPCC. However, to 
minimize fragmentation of the machine, jobs using up to 32 nodes were limited to a total ex- 
ecution time of 30 minutes. One possibility was to break up the run into 30 minute portions, 
but this introduces additional overhead and complications; however, new jobs start instantly 
because they need not be queued. A better option was to run jobs on a 64-node partition, 
task farming four 16-node jobs to run simultaneously. There was a 12 hour limit to 64-512 
node jobs (6 hours during the week), but there was often a long wait in the queues. If the 
efficient usage of billed time was a significant concern, sixteen 2-node jobs would complete 
in 8.0 hours. However, during the week this meant restarting halfway through and waiting 
in the queue again. Similar comments apply to the Hitachi SR2201, although its queues 
were limited to 8 hours maximum. The extra administrative overheads involved and the 
billed usage means that we usually concentrated computation on the serial workstations. 
However, parallel execution becomes more attractive with larger simulations. 

The parallel efficiency of DPD with 6400 particles is good only for a modest number 
of processor nodes. This is particularly true of the more modern parallel machines such 
as the Cray T3E and SGI Origin2000, which are proportionally faster in processing than 
communicating when compared with their older counterparts. Much better parallel efficiency 
has been observed with larger simulations. The Cray T3D shows an unusual increase in 
efficiency when going from a serial calculation to a 2-node parallel calculation with 6400 
particles; this could be explained by any number of hardware-specific arguments. We should 
note that the results for the Origin2000 include the effect of sharing the machine with other 
users, unlike all the other machines whose results appear in Tables |VIII| and [TX|, for which 
each node was dedicated to our calculations. 

We decided to write the main simulation program in C/C++ as opposed to Fortran. 
This choice was made because C/C++ were believed to be the most appropriate languages 
for dealing with DPD simulations which consist of a large amount of book-keeping wrapped 
around fairly simple computations. C/C++ and Fortran are highly portable to different 
computer architectures, and although well-written Fortran is more efficient on vector ma- 
chines, for almost all other situations they are of similar speed, given equally good compilers. 
The use of vector machines (such as the Hitachi SR2201) was not anticipated when this work 
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on DPD began several years ago. Furthermore, it was not believed that the basic algorithm 
would vectorize well, due to the short vector length in typical computations. Large programs 
are easier to maintain in C/C++ than in Fortran, although the increasingly well-supported 
Fortran 90 and 95 make the difference less significant. 

Finally, we comment on our findings in tuning the message passing interface (MPI) calls 
for the Cray T3D. In our simulations, it was found that blocking calls (sends and receives) 
were faster than non-blocking calls and were easier to use correctly. Furthermore, better 
scaling was achieved by sending the size of a variable-size message in a separate message 
rather than probing incoming messages to determine their size. Finally, using derived data 
types to remove unneeded data from messages was slower than sending everything. 
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TABLE I. Model parameters for the viscosity measurements of the lower-viscosity fluid. 
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TABLE II. Model parameters for the viscosity measurements of the higher-viscosity fluid. 
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TABLE III. Model parameters used for the differing-viscosity immiscible fluid simulations. 
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TABLE IV. Scaling exponents for two-dimensional differing-viscosity fluids as a function of 
viscous phase fraction, divided into "early" and "late time". 
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U.OD41 It u.uu / 


0.7 


0.85 ±0.06 to 2.62 ±0.10 


0.363 ±0.012 


0.8 


1.08 ±0.04 to 2.51 ±0.06 


0.406 ±0.017 


0.9 


1.76 ± 0.08 to 2.70 ± 0.03 


U.OZ zt u.uz 


TABLE V. Scaling exponents for three-dimensional differing-viscosity fluids 


dh d IUIlCTilOIl OI 


viscous 


phase fraction, divided into "early" and "late time". 






Approximate range of log±ot 


p 


0.2 


to 1.10 ±0.04 


0.123 ± 0.016 


0.25 


0.02 ±0.04 to 1.04 ±0.06 


c\ 1 o i o no 

0.13 ± 0.02 


0.3 


0.02 ±0.08 to 1.00 ±0.05 


oior - i omo 

0.135 ± 0.012 


0.4 


-0.015 ±0.010 to 0.85 ±0.03 


r\ ~i ~i a i n ni « 

0.114 ± 0.017 


0.5 


0.00 ±0.03 to 0.90 ±0.08 


o i o i o no 

0.13 ± 0.03 


0.05 


0.19 ±0.19 to 2.83 ±0.14 


0.283 ± 0.019 


0.1 


0.12 ±0.16 to 2.81 ±0.18 


0.304 ± 0.010 


0.15 


0.24 ±0.20 to 2.71 ±0.18 


304 ± 010 


0.2 


0.52 ±0.24 to 2.71 ±0.15 


0.337 ±0.012 


0.25 


1.0 ±0.3 to 2.76 ±0.13 


0.39 ± 0.04 


0.3 


0.71 ±0.15 to 2.53 ±0.22 


0.367 ±0.016 


0.4 


1.11 ±0.16 to 2.5 ±0.2 


0.415 ±0.019 


0.5 


1.4 ±0.3 to 2.50 ±0.09 


0.47 ± 0.04 



TABLE VI. Scaling exponents for two-dimensional equal-viscosity fluids as a function of mi- 
nority phase fraction, divided into "early" and "late time" . 
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9 


Approximate range of logiot 




P 


0.1 


to 1.41 ± 0.07 




0.028 ± 0.013 


0.2 


r\ r\-\ i r\ r\ a i r\ r\r\ i r\ r\ A 

0.01 ± 0.04 to 0.90 ± 0.04 




0.026 ± 0.007 


0.3 


0.00 ± 0.04 to 0.73 ± 0.05 




0.026 ± 0.007 


0.4 


to 0.61 ± 0.05 




0.023 ± 0.005 


0.5 


to 60 ± 05 




029 ± 006 


0.1 


1.35 ±0.11 to 2.48 ±0.18 




0.43 ± 0.08 


0.2 


0.52 ± 0.04 to 2.68 ± 0.02 




n q«9 -i- n nns 

u.uUi in u.uuo 


0.3 


0.43 ± 0.04 to 2.53 ± 0.04 




0.366 ± 0.009 


0.4 


0.38 ±0.02 to 2.47 ±0.03 




0.369 ± 0.006 


0.5 


0.32 ± 0.04 to 2.48 ± 0.02 




n oca J- n (1(17 


TABLE VII. Scaling exponents for three-dimensional equal-viscosity fluids 


as f\ f n n ct,\ on of 


minority phase fraction, divided into "early" and "late time" . 








Number 


Elapsed time 


Parallel 


Machine 


of nodes 


per node 


em cue n cy 


DEC 3000/400 AXP 


1 


624 




(133 MHz Alpha EV4) 








Linux PC 


1 


180 




(350 MHz Intel Pentium II) 








SGI Indigo2 


1 


200 




(195 MHz MIPS R10000) 








EPCC Cray T3D 


1 


1254 


1 an 

1.00 


(512 x 150 MHz Alpha EV4) 


2 


575 


1.09 




4 


354 


0.89 




8 


251 


0.62 




16 


206 


0.38 




32 


214 


0.18 


HPCF Hitachi SR2201 


1 


1202 


1.00 


(256 x 150 MHz HARP-IE) 


2 


634 


0.95 




4 


371 


0.81 




8 


255 


0.59 




16 


212 


0.35 




32 


243 


0.15 



TABLE VIII. Elapsed time (in seconds) per node and parallel efficiency of various computers 
for the first 1000 time steps of a 6400 particle simulation in two dimensions. These computers were 
used to calculate the results in this paper. 
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Number 


Elapsed time 


Parallel 


Machine 


or nodes 


per node 


efficiency 


bGl Octane 


1 


152 


1.00 


[ZOU IViriZ Mlro xtiUUUU J 


o 
z 


Q7 
O I 


n qq 
U.oo 


ubAK uray 1311j-120011j 


1 


1 A O 

14,5 


1.00 


{5(0 x bUU IViriz Alpha HvVoJ 


2 


9b 


0. /4 




4 


75 


0.48 




8 


67 


0.27 




16 


70 


0.13 


CSAR SGI Origin2000 


1 


133 


1.00 


(16 x 250 MHz MIPS R10000) 


2 


81 


0.83 




4 


69 


0.49 




8 


60 


0.28 



TABLE IX. Elapsed time (in seconds) per node and parallel efficiency of various computers for 
the first 1000 time steps of a 6400 particle simulation in two dimensions. These computers were 
not used to calculate the results in this paper. 
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FIGURES 



FIG. 1. Time-evolution of five sample simulations of two-dimensional differing- viscosity fluids, 
each simulation having a different viscous phase fraction (varying from (ft = 0.1 through (ft = 0.5). 
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FIG. 2. Time-evolution of five sample simulations of two-dimensional differing- viscosity fluids, 
each simulation having a different viscous phase fraction (varying from cf> = 0.5 through <p = 0.9). 
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FIG. 3. Time-evolution of five sample simulations of three-dimensional differing-viscosity fluids, 
each simulation having a different viscous phase fraction (varying from (f> = 0.1 through <j> = 0.5). 
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FIG. 4. Time-evolution of five sample simulations of three-dimensional differing-viscosity fluids, 
each simulation having a different viscous phase fraction (varying from cf> = 0.5 through <p = 0.9). 
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FIG. 5. Scaling exponents for two-dimensional differing-viscosity fluids as a function of viscous 
phase fraction. Circles indicate "early time" and horizontal marks "late time" ; error bars are 68% 
confidence intervals. 
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FIG. 6. Scaling exponents for three-dimensional differing-viscosity fluids as a function of viscous 
phase fraction. Circles indicate "early time" and horizontal marks "late time" ; error bars are 68% 
confidence intervals. 
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FIG. 7. Time-evolution of five sample simulations of two-dimensional equal- viscosity fluids, 
each simulation having a different minority phase fraction (varying from (f> = 0.1 through (p = 0.5). 
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FIG. 8. Time-evolution of five sample simulations of three-dimensional equal-viscosity fluids, 
each simulation having a different minority phase fraction (varying from (f> = 0.1 through (ft = 0.5). 
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FIG. 9. Scaling exponents for two-dimensional equal-viscosity fluids as a function of minority 
phase fraction. Circles indicate "early time" and horizontal marks "late time" ; error bars are 68% 
confidence intervals. 
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FIG. 10. Scaling exponents for three-dimensional equal- viscosity fluids as a function of minority 
phase fraction. Circles indicate "early time" and horizontal marks "late time" ; error bars are 68% 
confidence intervals. 
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FIG. 11. Differing-phase radial distribution function for a three-dimensional equal- viscosity 
fluid ((f) = 0.5) at t = 13.9. The unit of length for r is specific to the particular set of model 
parameters used. 
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